# load packages
library(expss)
library(knitr)
library(sandwich)
library(lmtest)
library(plyr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(arm)
library(foreign)
library(lattice)
library(MASS)
library(stats)
library(base)
library(car)
library(xtable)
library(broom)
library(stringr)
library(grid)
library(gridExtra)
library(lme4)
library(stringr)
library(plm)
library(estimatr)
library(lmtest)
library(multiwayvcov)
library(modelr)
library(readr)
library(haven)
library(stargazer)
library(brms)
library(jtools)
library(sjPlot)
library(sjlabelled)
library(sjmisc)
library(ggeffects)
library(clubSandwich)
library(texreg)
library(fixest)
library(cregg)
library(cjoint)
library(clusterSEs)
library(knitr)
library(kableExtra)
library(magrittr)
library(tidyverse)
library(rstatix)
library(ggpubr)
library(gtsummary)
library(gt)

rm(list = ls()) 

options(scipen=999)#turn off scientific notation or price of planes will be in scientific notation and the code for conjoint won't work.

multi <- read.csv("multi.csv") #load data


# group is the treatment variable. news is the reference level, then there are two other groups (generic info and specific conjoint info)
multi$group <- as.factor(multi$group)
multi$group <- relevel(multi$group, ref = "News")

multi$group <-factor(multi$group, levels = c("News",  "Specific Information", "Generic Information"))

# base models fairness is the DV 
mod1 <- lm(agree_1 ~ group, data=multi)
pred1 <- ggpredict(mod1, terms = c("group"))#predicted values

#this is figure 1
plot(pred1) + 
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high, x=x, width = 0.1))+
  labs(x = "", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "DV: The company's decision to introduce the new technology is fair") + 
  theme(text = element_text(size = 18)) + ylim(1,5) + coord_flip() 

# Here DV is whether you'd do the same thing if you were CEO
mod2 <- lm(agree_2 ~ group, data=multi)
pred2 <- ggpredict(mod2, terms = c("group"))
# figure A9
pdf(file = "figureA9.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(pred2) + 
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high, x=x, width = 0.1))+
  labs(x = "", y="Agree (1 - Strongly disagree to 5 - Strongly agree)",
       title= "DV: If you were CEO you'd make same choice")  + 
  theme(text = element_text(size = 18)) + ylim(1,5) + coord_flip()
dev.off()

# here the DVs are support for different policies in response to automation
mod3 <- lm(policy_socialspend ~ group, data=multi)
pred3 <- ggpredict(mod3, terms = c("group"))
#figure A1
pdf(file = "figureA1.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(pred3) + 
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high, x=x, width = 0.1))+
  labs(x = "", y="Agree (1 - Strongly disagree to 5 - Strongly agree)",
       title= "DV: Social Spending")  + 
  theme(text = element_text(size = 18)) + ylim(1,5) + coord_flip()
dev.off()

mod4 <- lm(policy_basicincome ~ group, data=multi)
pred4 <- ggpredict(mod4, terms = c("group"))
#figure A2
pdf(file = "figureA2.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(pred4) + 
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high, x=x, width = 0.1))+
  labs(x = "", y="Agree (1 - Strongly disagree to 5 - Strongly agree)",
       title= "DV: Basic Income") + 
  theme(text = element_text(size = 18)) + ylim(1,5) + coord_flip()
dev.off()

mod5 <- lm(policy_jobguarantee ~ group, data=multi)
#figure A3
pdf(file = "figureA3.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
pred5 <- ggpredict(mod5, terms = c("group"))
plot(pred5) + 
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high, x=x, width = 0.1))+
  labs(x = "", y="Agree (1 - Strongly disagree to 5 - Strongly agree)",
       title= "DV: Job Guarantee") + 
  theme(text = element_text(size = 18))+ ylim(1,5) + coord_flip()
dev.off()

mod6 <- lm(policy_unskilled ~ group, data=multi)
#figure A4
pdf(file = "figureA4.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
pred6 <- ggpredict(mod6, terms = c("group"))
plot(pred6) + geom_errorbar(aes(ymin=conf.low, ymax=conf.high, x=x, width = 0.1))+
  labs(x = "", y="Agree (1 - Strongly disagree to 5 - Strongly agree)",
       title= "DV: Unskilled") + 
  theme(text = element_text(size = 18)) + ylim(1,5)+ coord_flip()
dev.off()

mod7 <- lm(policy_skilled ~ group, data=multi)
pred7 <- ggpredict(mod7, terms = c("group"))
#figure A5
pdf(file = "figureA5.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(pred7) + 
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high, x=x, width = 0.1))+
  labs(x = "", y="Agree (1 - Strongly disagree to 5 - Strongly agree)",
       title= "DV: Skilled") + 
  theme(text = element_text(size = 18)) + ylim(1,5)+ coord_flip()
dev.off()

mod8 <- lm(policy_trade ~ group, data=multi)
pred8 <- ggpredict(mod8, terms = c("group"))
#figure A6
pdf(file = "figureA6.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(pred8) + 
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high, x=x, width = 0.1))+
  labs(x = "", y="Agree (1 - Strongly disagree to 5 - Strongly agree)",
       title= "DV: Trade") + 
  theme(text = element_text(size = 18)) + ylim(1,5)+ coord_flip()
dev.off()

mod9 <- lm(policy_reskill ~ group, data=multi)
pred9 <- ggpredict(mod9, terms = c("group"))
#figure A7
pdf(file = "figureA7.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(pred9) + 
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high, x=x, width = 0.1))+
  labs(x = "", y="Agree (1 - Strongly disagree to 5 - Strongly agree)",
       title= "DV: Reskill") + 
  theme(text = element_text(size = 18)) + ylim(1,5)+ coord_flip()
dev.off()

mod10 <- lm(policy_automationtax ~ group, data=multi)
pred10 <- ggpredict(mod10, terms = c("group"))
#figure A8
pdf(file = "figureA8.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(pred10) + 
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high, x=x, width = 0.1))+
  labs(x = "", y="Agree (1 - Strongly disagree to 5 - Strongly agree)",
       title= "DV: Tax automation") + 
  theme(text = element_text(size = 18)) + ylim(1,5)+ coord_flip()
dev.off()

## regression tables

#table A1
stargazer(mod1,mod2,
          type = "latex", column.labels=c("Fairness", "Would do same if CEO"),
          out="TableA1.tex",
          model.numbers = F, dep.var.labels  = c("",""),
          covariate.labels=c("Specific Information (ref. News)","Generic Information","Constant"), star.cutoffs = c(0.05, 0.01, 0.001),
          star.char=c("*", "**", "***"),
          omit.stat = c("f", "adj.rsq", "ser"),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"),notes.append = FALSE)

# table A2
stargazer(mod3,mod4,mod5,mod6,mod7,mod8,mod9,mod10,
          out="TableA2.tex",
          dep.var.labels  = c("","","","","","","",""),
          type = "latex", column.labels=c("Social Spending","Basic Income", "Job Guarantee", "Unskilled", "Skilled",
                                          "Trade", "Reskilling", "Tax Automation"),
          covariate.labels=c("Specific Information (ref. News)","Generic Information", "Constant"), star.cutoffs = c(0.05, 0.01, 0.001),
          star.char=c("*", "**", "***"),
          omit.stat = c("f", "adj.rsq", "ser"),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"),notes.append = FALSE)


##############
# conjoint
###############
#select only variables we'll use
multi2 <- dplyr::filter(multi, group=="Specific Information") #this is the group that saw the conjoint
multi2 <- dplyr::select(multi2, numlsafter, numlsafter1, numlsafter2, numlsafter3, numhsafter, numhsafter1, 
                        numhsafter2, numhsafter3, wagelsafter, wagelsafter1, wagelsafter2, wagelsafter3, 
                        wagehsafter, wagehsafter1, wagehsafter2, wagehsafter3, priceafter, priceafter1, 
                        priceafter2, priceafter3, product, product1, product2, product3, scen1_agree_1,
                        scen2_agree_1, scen3_agree_1, scen4_agree_1,scen1_agree_2,
                        scen2_agree_2, scen3_agree_2, scen4_agree_2, ResponseId, UserLanguage, country,
                        scen1_know1, scen1_know2, female, college, LR, self_know_automation_cat, know_automation_correct,
                        screen1)

#rename variables 
multi2 <- dplyr::rename(multi2, numlsafter_1 = numlsafter, numlsafter_2 = numlsafter1, numlsafter_3 = numlsafter2,
                        numlsafter_4 = numlsafter3, numhsafter_1 = numhsafter, numhsafter_2 = numhsafter1, 
                        numhsafter_3 = numhsafter2, numhsafter_4 = numhsafter3, wagelsafter_1 = wagelsafter,
                        wagelsafter_2 = wagelsafter1, wagelsafter_3 = wagelsafter2, wagelsafter_4 = wagelsafter3, 
                        wagehsafter_1 = wagehsafter, wagehsafter_2 =wagehsafter1, wagehsafter_3 = wagehsafter2,
                        wagehsafter_4 = wagehsafter3, priceafter_1 = priceafter, priceafter_2 = priceafter1, 
                        priceafter_3 = priceafter2, priceafter_4 = priceafter3, product_1 = product, 
                        product_2 = product1, product_3 = product2, product_4 = product3, 
                        scenagree_1_1 = scen1_agree_1,
                        scenagree_2_1 = scen2_agree_1,  scenagree_3_1 = scen3_agree_1, 
                        scenagree_4_1 = scen4_agree_1, scenagree_1_2 = scen1_agree_2,
                        scenagree_2_2 = scen2_agree_2, scenagree_3_2 = scen3_agree_2, 
                        scenagree_4_2 = scen4_agree_2, scenknow_1_1 = scen1_know1 , scenknow_1_2 = scen1_know2)

#reshape dataset so we can analyze it
multi2_reshape <- multi2 %>% 
  #first step is gathering all the variables that we want to use as our attributes;
  #because we have named each task iteration in this format, using contains will get each conjoint-related variable
  gather(key = task, value = score, contains("_1"), contains("_2"), contains("_3"), contains("_4")) %>% 
  #we need to separate these variables out by characteristics; this can be done using `separate`, because variable were meaningfully named using underscores 
  separate(task, c("variable", "iteration", "index")) %>% 
  #creating an indicator of variables that were a part of an index
  dplyr::mutate(index = ifelse(is.na(index), "", index))  %>%
  #create a variable specific column
  unite("variable_index", c("variable", "index"), sep = "") %>% 
  #spreading out the variables so that each meaningful variable (or attribute) is a column
  tidyr::spread(key = variable_index, value = score) %>% 
  #lastly, because the data are a mix of upper and lower case, I use clean_names to make the variables names follow a uniform convention (this is optional, but recommended)
  janitor::clean_names() 

#create factor variables

multi2_reshape <- multi2_reshape %>%
  mutate(female = case_when(
    female== 0 ~"Male",
    female == 1 ~ "Female"
  ))

multi2_reshape <- multi2_reshape %>%
  mutate(college = case_when(
    college== 0 ~"No College",
    college == 1 ~ "College"
  ))

#makes sure all of our variables of interest are factors
multi2_reshape$female <- as.factor(multi2_reshape$female)
multi2_reshape$college <- as.factor(multi2_reshape$college)
multi2_reshape$know_automation_correct <- dplyr::recode(multi2_reshape$know_automation_correct,`0`="Low",`1`= "High")


multi2_reshape$scenagree1 <- as.numeric(as.character(multi2_reshape$scenagree1))
multi2_reshape$scenagree2 <- as.numeric(as.character(multi2_reshape$scenagree2))
multi2_reshape$numhsafter <- as.factor(as.character(multi2_reshape$numhsafter))
multi2_reshape$numlsafter <- as.factor(as.character(multi2_reshape$numlsafter))
multi2_reshape$priceafter <- as.factor(as.character(multi2_reshape$priceafter))
multi2_reshape$product <- as.factor(as.character(multi2_reshape$product))
multi2_reshape$wagehsafter <- as.factor(as.character(multi2_reshape$wagehsafter))
multi2_reshape$wagelsafter <- as.factor(as.character(multi2_reshape$wagelsafter))
multi2_reshape$self_know_automation_cat <- as.factor(as.character(multi2_reshape$self_know_automation_cat))

#recode attributes
multi2_reshape$priceafter <- ifelse(multi2_reshape$priceafter=="100000000" | 
                                      multi2_reshape$priceafter== "25" | 
                                      multi2_reshape$priceafter== "25000"| 
                                      multi2_reshape$priceafter== "600", "Same price", 
                                    ifelse(multi2_reshape$priceafter=="80000000" | 
                                             multi2_reshape$priceafter== "20" | 
                                             multi2_reshape$priceafter== "20000"| 
                                             multi2_reshape$priceafter== "480", "20% cheaper", 
                                           ifelse(multi2_reshape$priceafter=="50000000" | 
                                                    multi2_reshape$priceafter== "12.5" | 
                                                    multi2_reshape$priceafter== "12500"| 
                                                    multi2_reshape$priceafter== "300", "50% cheaper",  NA)))

multi2_reshape$product <- ifelse(multi2_reshape$product=="Vaccine" | 
                                   multi2_reshape$product=="Vaccin", "Vaccine",
                                 ifelse(multi2_reshape$product=="smartphone" | 
                                          multi2_reshape$product=="Téléphone intelligent", "Smartphone",
                                        ifelse(multi2_reshape$product=="Automobile" | 
                                                 multi2_reshape$product=="Car", "Car",
                                               ifelse(multi2_reshape$product=="Avion" | 
                                                        multi2_reshape$product=="Plane", "Plane",  NA))))

multi2_reshape$wagehsafter <- ifelse(multi2_reshape$wagehsafter=="£ 75,000" | 
                                       multi2_reshape$wagehsafter=="$125,000", "$125,000",
                                     ifelse(multi2_reshape$wagehsafter=="£ 90,000" | 
                                              multi2_reshape$wagehsafter=="$150,000", "$150,000",NA))

multi2_reshape$wagelsafter <- ifelse(multi2_reshape$wagelsafter=="£ 13,000" | 
                                       multi2_reshape$wagelsafter=="$20,000"| 
                                       multi2_reshape$wagelsafter=="$30,000",
                                     "$20,000",
                                     ifelse(multi2_reshape$wagelsafter=="£ 17,000" | 
                                              multi2_reshape$wagelsafter=="$25,000"| 
                                              multi2_reshape$wagelsafter=="$40,000",
                                            "$25,000",NA))

#recode knowledge
multi2_reshape$scenknow1 <- ifelse(multi2_reshape$scenknow1=="400 before and ${e://Field/numbertotal2} after",
                                   "1","0")
multi2_reshape$scenknow2 <- ifelse(multi2_reshape$scenknow2=="${e://Field/price}%",
                                   "1","0")

multi2_reshape$scenknow_all <- ifelse(multi2_reshape$scenknow1=="1" & multi2_reshape$scenknow2=="1",
                                      "High Knowledge",
                                      ifelse((multi2_reshape$scenknow1=="1" & multi2_reshape$scenknow2=="0") |
                                               (multi2_reshape$scenknow2=="1" & multi2_reshape$scenknow1=="0"), "Medium Knowledge",
                                             "Low Knowledge"))

multi2_reshape$priceafter <- as.factor(as.character(multi2_reshape$priceafter))
multi2_reshape$priceafter <- relevel(multi2_reshape$priceafter, ref = c("Same price"))

#recode into factors
multi2_reshape$numhsafter <- as.factor(as.character(multi2_reshape$numhsafter))
multi2_reshape$numlsafter <- as.factor(as.character(multi2_reshape$numlsafter))
multi2_reshape$product <- as.factor(as.character(multi2_reshape$product))
multi2_reshape$wagehsafter <- as.factor(as.character(multi2_reshape$wagehsafter))
multi2_reshape$wagelsafter <- as.factor(as.character(multi2_reshape$wagelsafter))
multi2_reshape$screen1 <- as.factor(as.character(multi2_reshape$screen1))
multi2_reshape$scenknow_all <- as.factor(as.character(multi2_reshape$scenknow_all))
multi2_reshape$scenknow_all <- relevel(multi2_reshape$scenknow_all, ref = c("Low Knowledge"))
multi2_reshape$self_know_automation_cat <- relevel(multi2_reshape$self_know_automation_cat, ref = c("Low Self-Know"))
#labels
attr(multi2_reshape$numhsafter, "label") <- "Number of High Skilled - After"
attr(multi2_reshape$numlsafter, "label") <- "Number of Low Skilled - After"
attr(multi2_reshape$wagelsafter, "label") <- "Wage of Low Skilled - After"
attr(multi2_reshape$wagehsafter, "label") <- "Wage of High Skilled - After"
attr(multi2_reshape$priceafter, "label") <- "Price - After"
attr(multi2_reshape$product, "label") <- "Product"


#### models #####

#marginal means for DV fairness
#figure 3
plot(cregg::mm(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id), size=3) + theme(text = element_text(size = 15))   

mm <- cregg::mm(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
            priceafter + product, id= ~ response_id)

mm$outcome <- "fairness"

#table of marginal means table A4
mm %>%
  kable(format = "latex", booktabs = TRUE, label="appendix_fairness") %>%
  kable_classic(full_width = F, font_size=7, latex_options = c("hold_position")) %>%
  save_kable("TableA4.tex")

# amce for fairness
# figure 4
plot(cj(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id),
        size=3) + theme(text = element_text(size = 15))   
cj <- cj(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
     priceafter + product, id= ~ response_id)

cj$outcome <- "fairness"

## table A5

cj %>%
  kable(format = "latex", booktabs = TRUE, label="appendix_fairness_amce") %>%
  kable_classic(full_width = F, font_size=7, latex_options = c("hold_position")) %>%
  save_kable("TableA5.tex")

#marginal means for CEO
#figure A10
pdf(file = "figureA10.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(cregg::mm(multi2_reshape, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id), size=3) + theme(text = element_text(size = 15))  
dev.off()

# amce for CEO
#figure A11
pdf(file = "figureA11.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(cj(multi2_reshape, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id), size=3) + theme(text = element_text(size = 15)) 
dev.off()

#### Subgroup analysis ######
## fairness ###
## by objective knowledge

#marginal means
mm_by <- cregg::cj(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product + college, id= ~ response_id, estimate = "mm", by = ~scenknow_all)
plot(mm_by, group = "scenknow_all", vline = 0.0, size=3) + theme(text = element_text(size = 15)) 

#amces
amces <- cregg::cj(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product + college, id= ~ response_id, estimate = "amce", 
            by = ~scenknow_all)
#differences in amces
diff_amces <- cregg::cj(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product + college, id= ~ response_id,estimate = "amce_diff", 
                 by = ~scenknow_all)

#figure A14
pdf(file = "figureA14.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()

# only show high and low knowledge for main text figure 5
diff_amces <- filter(diff_amces, BY=="High Knowledge - Low Knowledge")
amces <- filter(amces, BY!="Medium Knowledge")
amces$BY <- droplevels(amces$BY)
levels(amces$BY) <- c("High Know", "Low Know")
diff_amces$BY[diff_amces$BY == "High Knowledge - Low Knowledge"] <- "High - Low" 

plot(rbind(amces,diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 18)) #figure 5

# table A6
knowledge_amces_low <- filter(amces, BY=="Low Know")
knowledge_amces_low <- dplyr::select(knowledge_amces_low,-scenknow_all)
knowledge_amces_low$outcome <- "fairness"
knowledge_amces_low$BY <- "Low Knowledge"

knowledge_amces_high <- filter(amces, BY=="High Know")
knowledge_amces_high <- dplyr::select(knowledge_amces_high,-scenknow_all)
knowledge_amces_high$outcome <- "fairness"
knowledge_amces_high$BY <- "High Knowledge"

knowledge_diff_amces <- diff_amces
knowledge_diff_amces <- dplyr::select(knowledge_diff_amces,-scenknow_all)
knowledge_diff_amces$outcome <- "fairness"

#table for low knowledge A6
knowledge_amces_low %>%
  kable(format = "latex", booktabs = TRUE, label="appendix_knowledge_amce_low") %>%
  kable_classic(full_width = F, font_size=7, latex_options = c("hold_position")) %>%
  save_kable("TableA6.tex")

#table for high knowledge A7
knowledge_amces_high %>%
  kable(format = "latex", booktabs = TRUE, label="appendix_knowledge_amce_high") %>%
  kable_classic(full_width = F, font_size=7, latex_options = c("hold_position")) %>%
  save_kable("TableA7.tex")

# table for difference between high and low knowledge A8
knowledge_diff_amces %>%
  kable(format = "latex", booktabs = TRUE, label="appendix_knowledge_diff_amce") %>%
  kable_classic(full_width = F, font_size=7, latex_options = c("hold_position")) %>%
  save_kable("TableA8.tex")

# calculate conditional MMs
mms <- cj(na.omit(multi2_reshape), scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
            priceafter + product + college, id= ~ response_id,  estimate = "mm", by = ~scenknow_all)
diff_mms <- cj(na.omit(multi2_reshape), scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                 priceafter + product + college, id= ~ response_id, estimate = "mm_diff", 
               by = ~scenknow_all)
plot(rbind(mms, diff_mms)) + ggplot2::facet_wrap(~BY, ncol = 3L)


## subjective self-knowledge
# fairness

#amces 
amces <- cj(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product, id= ~ response_id, estimate = "amce", 
            by = ~self_know_automation_cat)
#differences in amces
diff_amces <- cj(multi2_reshape, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product, id= ~ response_id,estimate = "amce_diff", 
                 by = ~self_know_automation_cat)

# figure A15
pdf(file = "figureA15.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()

# high and low self-knowledge only for main text plot figure 6
diff_amces <- filter(diff_amces, BY=="High Self-Know - Low Self-Know")
amces <- filter(amces, BY!="Medium Self-Know")
amces$BY <- droplevels(amces$BY)
levels(amces$BY) <- c("High Self-Know", "Low Self-Know")
diff_amces$BY[diff_amces$BY == "High Self-Know - Low Self-Know"] <- "High - Low" 

plot(rbind(amces,diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 18)) #figure 6

# tables 
self_knowledge_amces_low <- filter(amces, BY=="Low Self-Know")
self_knowledge_amces_low <- dplyr::select(self_knowledge_amces_low,-self_know_automation_cat)
self_knowledge_amces_low$outcome <- "fairness"
self_knowledge_amces_low$BY <- "Low Self-Knowledge"

self_knowledge_amces_high <- filter(amces, BY=="High Self-Know")
self_knowledge_amces_high <- dplyr::select(self_knowledge_amces_high,-self_know_automation_cat)
self_knowledge_amces_high$outcome <- "fairness"
self_knowledge_amces_high$BY <- "High Self-Knowledge"

self_knowledge_diff_amces <- diff_amces
self_knowledge_diff_amces <- dplyr::select(self_knowledge_diff_amces,-self_know_automation_cat)
self_knowledge_diff_amces$outcome <- "fairness"

#table for low self-knowledge A9
self_knowledge_amces_low %>%
  kable(format = "latex", booktabs = TRUE, label="appendix_self_knowledge_amce_low") %>%
  kable_classic(full_width = F, font_size=7, latex_options = c("hold_position")) %>%
  save_kable("TableA9.tex")

#table for high self-knowledge A10
self_knowledge_amces_high %>%
  kable(format = "latex", booktabs = TRUE, label="appendix_self_knowledge_amce_high") %>%
  kable_classic(full_width = F, font_size=7, latex_options = c("hold_position")) %>%
  save_kable("TableA10.tex")

# table for difference between high and low A11
self_knowledge_diff_amces %>%
  kable(format = "latex", booktabs = TRUE, label="appendix_self_knowledge_diff_amce") %>%
  kable_classic(full_width = F, font_size=7, latex_options = c("hold_position")) %>%
  save_kable("TableA11.tex")


### CEO

## CEO ###
## by objective knowledge
#marginal means
mm_by <- cj(multi2_reshape, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product + college, id= ~ response_id, estimate = "mm", by = ~scenknow_all)
plot(mm_by, group = "scenknow_all", vline = 0.0, size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
#amces
amces <- cj(multi2_reshape, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product + college, id= ~ response_id, estimate = "amce", 
            by = ~scenknow_all)
#difference in amces
diff_amces <- cj(multi2_reshape, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product + college, id= ~ response_id,estimate = "amce_diff", 
                 by = ~scenknow_all)

# figure A12
pdf(file = "figureA12.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()

# calculate conditional MMs
mms <- cj(multi2_reshape, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
            priceafter + product, id= ~ response_id,  estimate = "mm", by = ~scenknow_all)
diff_mms <- cj(multi2_reshape, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                 priceafter + product, id= ~ response_id, estimate = "mm_diff", 
               by = ~scenknow_all)
plot(rbind(mms, diff_mms), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 

## self-knowledge
#amces
amces <- cj(multi2_reshape, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product, id= ~ response_id, estimate = "amce", 
            by = ~self_know_automation_cat)
#difference in amces
diff_amces <- cj(multi2_reshape, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product, id= ~ response_id,estimate = "amce_diff", 
                 by = ~self_know_automation_cat)
# figure A13
pdf(file = "figureA13.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()


#### compare best conjoint scenario to worst conjoint scenario

#fairness

# calculate interaction MMs
mms <- cj(multi2_reshape, scenagree1 ~ numhsafter + wagehsafter + wagelsafter +
            product, id= ~ response_id, estimate = "mm", by = ~ priceafter + numlsafter)
plot(mms, vline = 0.5) + ggplot2::facet_wrap(~BY, ncol = 3L)

mms

best_worst <- dplyr::filter(mms, (BY=="50% cheaper***150" & level == "350") | (BY=="Same price***50" & level == "200"))

best_worst$outcome <- "fairness"

best_worst <- dplyr::mutate(best_worst, scen = ifelse(BY=="Same price***50", "Worst-case", "Best-case"))



# table A3
best_worst <- dplyr::select(best_worst, -BY)
best_worst %>%
  kable(format = "latex", booktabs = TRUE, label="best_worst_mm") %>%
  kable_classic(full_width = F, font_size=7, latex_options = c("hold_position")) %>%
  save_kable("TableA3.tex")


## 50% cheaper,150 and hs 350 3.6 (3.53, 3.66)
## same price, 50 and hs 200 3.09 (3.02, 3.16)

## by country: tables A12 through A48

### by country
# fairness
mod1_US <- lm(agree_1 ~ group, data = subset(multi, 
                                             country == "US"))
pred1_US <- ggpredict(mod1_US, terms = c("group"))
plot_US <- plot(pred1_US) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "US")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)

mod1_AUS <- lm(agree_1 ~ group, data = subset(multi, 
                                              country == "Australia"))
pred1_AUS <- ggpredict(mod1_AUS, terms = c("group"))
plot_AUS <- plot(pred1_AUS) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Australia")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)


mod1_CAN <- lm(agree_1 ~ group, data = subset(multi, 
                                              country == "Canada"))
pred1_CAN <- ggpredict(mod1_CAN, terms = c("group"))
plot_CAN <- plot(pred1_CAN) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Canada")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)


mod1_UK <- lm(agree_1 ~ group, data = subset(multi, 
                                             country == "UK"))
pred1_UK <- ggpredict(mod1_UK, terms = c("group"))
plot_UK <- plot(pred1_UK) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "UK")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)

grid.arrange(plot_US, plot_UK, plot_CAN,plot_AUS, ncol=2)

## fairness tables

stargazer(mod1_AUS,mod1_CAN,mod1_UK,mod1_US,
          type = "latex", 
          out="TableA12.tex",dep.var.labels=c("Fairness"),column.labels=c("Australia", "Canada", "UK", "US"),
          covariate.labels=c("Specific Information (ref. News)","Generic Information", "Constant"), star.cutoffs = c(0.05, 0.01, 0.001),
          star.char=c("*", "**", "***"),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"),omit.stat = c("f", "adj.rsq", "ser"),notes.append = FALSE)

## CEO

mod2_US <- lm(agree_2 ~ group, data = subset(multi, 
                                             country == "US"))
pred2_US <- ggpredict(mod2_US, terms = c("group"))
plot2_US <- plot(pred2_US) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "US")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)

mod2_AUS <- lm(agree_2 ~ group, data = subset(multi, 
                                              country == "Australia"))
pred2_AUS <- ggpredict(mod2_AUS, terms = c("group"))
plot2_AUS <- plot(pred2_AUS) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Australia")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)


mod2_CAN <- lm(agree_2 ~ group, data = subset(multi, 
                                              country == "Canada"))
pred2_CAN <- ggpredict(mod2_CAN, terms = c("group"))
plot2_CAN <- plot(pred2_CAN) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Canada")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)


mod2_UK <- lm(agree_2 ~ group, data = subset(multi, 
                                             country == "UK"))
pred2_UK <- ggpredict(mod2_UK, terms = c("group"))
plot2_UK <- plot(pred2_UK) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "UK")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)


grid.arrange(plot_US, plot_UK, plot_CAN,plot_AUS, ncol=2)

stargazer(mod2_AUS,mod2_CAN,mod2_UK,mod2_US,
          out="TableA13.tex",
          type = "latex", dep.var.labels=c("Would do the same if CEO"),column.labels=c("Australia", "Canada", "UK", "US"),
          covariate.labels=c("Specific Information (ref. News)","Generic Information", "Constant"), star.cutoffs = c(0.05, 0.01, 0.001),
          star.char=c("*", "**", "***"),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"),omit.stat = c("f", "adj.rsq", "ser"),notes.append = FALSE)

##  policies
# social spending
m3US <- lm(policy_socialspend ~ group, data = subset(multi, 
                                                           country == "US"))
pred3US <- ggpredict(m3US, terms = c("group"))
plot_US <- plot(pred3US) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "US")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)

m3au <- lm(policy_socialspend ~ group, data = subset(multi, 
                                                            country == "Australia"))
pred3AUS <- ggpredict(m3au, terms = c("group"))
plot_AUS <- plot(pred3AUS) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Australia")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)


m3CAN <- lm(policy_socialspend ~ group, data = subset(multi, 
                                                            country == "Canada"))
pred3CAN <- ggpredict(m3CAN, terms = c("group"))
plot_CAN <- plot(pred3CAN) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Canada")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)


m3UK <- lm(policy_socialspend ~ group, data = subset(multi, 
                                                           country == "UK"))
pred3UK <- ggpredict(m3UK, terms = c("group"))
plot_UK <- plot(pred3UK) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "UK")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)

# policy_basicincome
m4US <- lm(policy_basicincome ~ group, data = subset(multi, 
                                                           country == "US"))
pred4US <- ggpredict(m4US, terms = c("group"))
plot_US <- plot(pred4US) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "US")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)

m4au <- lm(policy_basicincome ~ group, data = subset(multi, 
                                                            country == "Australia"))
pred4AUS <- ggpredict(m4au, terms = c("group"))
plot_AUS <- plot(pred4AUS) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Australia")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)


m4CAN <- lm(policy_basicincome ~ group, data = subset(multi, 
                                                            country == "Canada"))
pred4CAN <- ggpredict(m4CAN, terms = c("group"))
plot_CAN <- plot(pred4CAN) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Canada")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)


m4UK <- lm(policy_basicincome ~ group, data = subset(multi, 
                                                           country == "UK"))
pred4UK <- ggpredict(m4UK, terms = c("group"))
plot_UK <- plot(pred4UK) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "UK")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)

# policy_jobguarantee

m5US <- lm(policy_jobguarantee ~ group, data = subset(multi, 
                                                            country == "US"))
pred5US <- ggpredict(m5US, terms = c("group"))
plot_US <- plot(pred5US) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "US")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)

m5au <- lm(policy_jobguarantee ~ group, data = subset(multi, 
                                                             country == "Australia"))
pred5AUS <- ggpredict(m5au, terms = c("group"))
plot_AUS <- plot(pred5AUS) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Australia")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)


m5CAN <- lm(policy_jobguarantee ~ group, data = subset(multi, 
                                                             country == "Canada"))
pred5CAN <- ggpredict(m5CAN, terms = c("group"))
plot_CAN <- plot(pred5CAN) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Canada")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)


m5UK <- lm(policy_jobguarantee ~ group, data = subset(multi, 
                                                            country == "UK"))
pred5UK <- ggpredict(m5UK, terms = c("group"))
plot_UK <- plot(pred5UK) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "UK")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)

# policy_unskilled
m6US <- lm(policy_unskilled ~ group, data = subset(multi, 
                                                         country == "US"))
pred6US <- ggpredict(m6US, terms = c("group"))
plot_US <- plot(pred6US) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "US")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)

m6au <- lm(policy_unskilled ~ group, data = subset(multi, 
                                                          country == "Australia"))
pred6AUS <- ggpredict(m6au, terms = c("group"))
plot_AUS <- plot(pred6AUS) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Australia")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)


m6CAN <- lm(policy_unskilled ~ group, data = subset(multi, 
                                                          country == "Canada"))
pred6CAN <- ggpredict(m6CAN, terms = c("group"))
plot_CAN <- plot(pred6CAN) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Canada")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)


m6UK <- lm(policy_unskilled ~ group, data = subset(multi, 
                                                         country == "UK"))
pred6UK <- ggpredict(m6UK, terms = c("group"))
plot_UK <- plot(pred6UK) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "UK")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)

# policy_skilled
m7US <- lm(policy_skilled ~ group, data = subset(multi, 
                                                       country == "US"))
pred7US <- ggpredict(m7US, terms = c("group"))
plot_US <- plot(pred7US) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "US")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)

m7au <- lm(policy_skilled ~ group, data = subset(multi, 
                                                        country == "Australia"))
pred7AUS <- ggpredict(m7au, terms = c("group"))
plot_AUS <- plot(pred7AUS) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Australia")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)


m7CAN <- lm(policy_skilled ~ group, data = subset(multi, 
                                                        country == "Canada"))
pred7CAN <- ggpredict(m7CAN, terms = c("group"))
plot_CAN <- plot(pred7CAN) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Canada")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)


m7UK <- lm(policy_skilled ~ group, data = subset(multi, 
                                                       country == "UK"))
pred7UK <- ggpredict(m7UK, terms = c("group"))
plot_UK <- plot(pred7UK) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "UK")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)

# policy_trade
m8US <- lm(policy_trade ~ group, data = subset(multi, 
                                                     country == "US"))
pred8US <- ggpredict(m8US, terms = c("group"))
plot_US <- plot(pred8US) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "US")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)

m8au <- lm(policy_trade ~ group, data = subset(multi, 
                                                      country == "Australia"))
pred8AUS <- ggpredict(m8au, terms = c("group"))
plot_AUS <- plot(pred8AUS) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Australia")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)


m8CAN <- lm(policy_trade ~ group, data = subset(multi, 
                                                      country == "Canada"))
pred8CAN <- ggpredict(m8CAN, terms = c("group"))
plot_CAN <- plot(pred8CAN) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Canada")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)


m8UK <- lm(policy_trade ~ group, data = subset(multi, 
                                                     country == "UK"))
pred8UK <- ggpredict(m8UK, terms = c("group"))
plot_UK <- plot(pred8UK) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "UK")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)


# policy_reskill
m9US <- lm(policy_reskill ~ group, data = subset(multi, 
                                                       country == "US"))
pred9US <- ggpredict(m9US, terms = c("group"))
plot_US <- plot(pred9US) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "US")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)

m9au <- lm(policy_reskill ~ group, data = subset(multi, 
                                                        country == "Australia"))
pred9AUS <- ggpredict(m9au, terms = c("group"))
plot_AUS <- plot(pred9AUS) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Australia")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)


m9CAN <- lm(policy_reskill ~ group, data = subset(multi, 
                                                        country == "Canada"))
pred9CAN <- ggpredict(m9CAN, terms = c("group"))
plot_CAN <- plot(pred9CAN) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Canada")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)


m9UK <- lm(policy_reskill ~ group, data = subset(multi, 
                                                       country == "UK"))
pred9UK <- ggpredict(m9UK, terms = c("group"))
plot_UK <- plot(pred9UK) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "UK")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)


# policy_automationtax
m10US <- lm(policy_automationtax ~ group, data = subset(multi, 
                                                              country == "US"))
pred10US <- ggpredict(m10US, terms = c("group"))
plot_US <- plot(pred10US) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "US")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)

m10au <- lm(policy_automationtax ~ group, data = subset(multi, 
                                                               country == "Australia"))
pred10AUS <- ggpredict(m10au, terms = c("group"))
plot_AUS <- plot(pred10AUS) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Australia")+
  theme(text = element_text(size = 12), legend.position = "none") +  ylim(2.5,4.4)


m10CAN <- lm(policy_automationtax ~ group, data = subset(multi, 
                                                               country == "Canada"))
pred10CAN <- ggpredict(m10CAN, terms = c("group"))
plot_CAN <- plot(pred10CAN) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "Canada")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)


m10UK <- lm(policy_automationtax ~ group, data = subset(multi, 
                                                              country == "UK"))
pred10UK <- ggpredict(m10UK, terms = c("group"))
plot_UK <- plot(pred10UK) + 
  labs(x = "Group", y="Agree (1 - Strongly disagree to 5 - Strongly agree)", title= "UK")+
  theme(text = element_text(size = 12), legend.position = "none") + ylim(2.5,4.4)

##Australia
stargazer(m3au,m4au,m5au,m6au,m7au,m8au,m9au,m10au,dep.var.labels.include = F,
          type = "latex", 
          out="TableA14.tex",column.labels=c("Social Spending","Basic Income", "Job Guarantee", "Unskilled", "Skilled",
                                          "Trade", "Reskilling", "Tax Automation"),
          covariate.labels=c("Specific Information (ref. News)","Generic Information", "Constant"), star.cutoffs = c(0.05, 0.01, 0.001),
          star.char=c("*", "**", "***"),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"),omit.stat = c("f", "adj.rsq", "ser"),notes.append = FALSE)

##Canada
stargazer(m3CAN,m4CAN,m5CAN,m6CAN,m7CAN,m8CAN,m9CAN,m10CAN,dep.var.labels.include = F,
          type = "latex", 
          out="TableA15.tex",column.labels=c("Social Spending","Basic Income", "Job Guarantee", "Unskilled", "Skilled",
                                          "Trade", "Reskilling", "Tax Automation"),
          covariate.labels=c("Specific Information (ref. News)","Generic Information", "Constant"), star.cutoffs = c(0.05, 0.01, 0.001),
          star.char=c("*", "**", "***"),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"),omit.stat = c("f", "adj.rsq", "ser"),notes.append = FALSE)

##UK
stargazer(m3UK,m4UK,m5UK,m6UK,m7UK,m8UK,m9UK,m10UK,dep.var.labels.include = F,
          type = "latex", 
          out="TableA16.tex",column.labels=c("Social Spending","Basic Income", "Job Guarantee", "Unskilled", "Skilled",
                                          "Trade", "Reskilling", "Tax Automation"),
          covariate.labels=c("Specific Information (ref. News)","Generic Information", "Constant"), star.cutoffs = c(0.05, 0.01, 0.001),
          star.char=c("*", "**", "***"),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"),omit.stat = c("f", "adj.rsq", "ser"),notes.append = FALSE)

##US
stargazer(m3US,m4US,m5US,m6US,m7US,m8US,m9US,m10US,dep.var.labels.include = F,
          type = "latex", 
          out="TableA17.tex",column.labels=c("Social Spending","Basic Income", "Job Guarantee", "Unskilled", "Skilled",
                                          "Trade", "Reskilling", "Tax Automation"),
          covariate.labels=c("Specific Information (ref. News)","Generic Information", "Constant"), star.cutoffs = c(0.05, 0.01, 0.001),
          star.char=c("*", "**", "***"),
          notes = c("* p<0.05; ** p<0.01; *** p<0.001"),omit.stat = c("f", "adj.rsq", "ser"),notes.append = FALSE)


## plot of policies

## US
summ <- multi %>%
  dplyr::filter(country=="US") %>%
  dplyr::group_by(group) %>%
  dplyr::summarize(socialsupport=mean(policy_socialspend, na.rm=T),
                   basicincome=mean(policy_basicincome, na.rm=T),
                   jobguarantee= mean(policy_jobguarantee, na.rm=T),
                   unskilled=mean(policy_unskilled, na.rm=T),
                   skilled=mean(policy_skilled, na.rm=T),
                   trade=mean(policy_trade, na.rm=T),
                   reskill=mean(policy_reskill, na.rm=T),
                   tax=mean(policy_automationtax, na.rm=T))

summ_long <- tidyr::gather(summ, policy_option, support, socialsupport: tax, 
                           factor_key=T)

US <- ggplot(summ_long, 
             aes(x=factor(policy_option), y=support, fill=group)) + 
  geom_bar(position="dodge", stat="identity") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 60, hjust=1, size=13),
        axis.text.y = element_text(size=14),
        axis.title.x = element_text(size=14),
        axis.title.y = element_text(size=14),
        plot.title = element_text(size=14, hjust=0.5, face="bold"), legend.position = "none") +
  ggtitle("US") +
  ylim(0,5) +
  xlab("") +
  ylab("Mean Support") +
  scale_x_discrete(labels= c("Social support", "Basic Income", "Job Guarantee", "Unskilled", "Skilled",
                             "Trade", "Reskill", "Tax automation"))

#australia

summ <- multi%>%
  dplyr::filter(country=="Australia") %>%
  dplyr::group_by(group) %>%
  dplyr::summarize(socialsupport=mean(policy_socialspend, na.rm=T),
                   basicincome=mean(policy_basicincome, na.rm=T),
                   jobguarantee= mean(policy_jobguarantee, na.rm=T),
                   unskilled=mean(policy_unskilled, na.rm=T),
                   skilled=mean(policy_skilled, na.rm=T),
                   trade=mean(policy_trade, na.rm=T),
                   reskill=mean(policy_reskill, na.rm=T),
                   tax=mean(policy_automationtax, na.rm=T))

summ_long <- tidyr::gather(summ, policy_option, support, socialsupport: tax, 
                           factor_key=T)

AUS <- ggplot(summ_long, 
              aes(x=factor(policy_option), y=support, fill=group)) + 
  geom_bar(position="dodge", stat="identity") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 60, hjust=1, size=13),
        axis.text.y = element_text(size=14),
        axis.title.x = element_text(size=14),
        axis.title.y = element_text(size=14),
        plot.title = element_text(size=14, hjust=0.5, face="bold"), legend.position = c(0.5,0.9), 
        legend.direction="horizontal") +
  ggtitle("Australia") +
  ylim(0,5) +
  xlab("") +
  ylab("Mean Support") +
  scale_x_discrete(labels= c("Social support", "Basic Income", "Job Guarantee", "Unskilled", "Skilled",
                             "Trade", "Reskill", "Tax automation"))


## UK
summ <- multi %>%
  dplyr::filter(country=="UK") %>%
  dplyr::group_by(group) %>%
  dplyr::summarize(socialsupport=mean(policy_socialspend, na.rm=T),
                   basicincome=mean(policy_basicincome, na.rm=T),
                   jobguarantee= mean(policy_jobguarantee, na.rm=T),
                   unskilled=mean(policy_unskilled, na.rm=T),
                   skilled=mean(policy_skilled, na.rm=T),
                   trade=mean(policy_trade, na.rm=T),
                   reskill=mean(policy_reskill, na.rm=T),
                   tax=mean(policy_automationtax, na.rm=T))

summ_long <- tidyr::gather(summ, policy_option, support, socialsupport: tax, 
                           factor_key=T)

UK <- ggplot(summ_long, 
             aes(x=factor(policy_option), y=support, fill=group)) + 
  geom_bar(position="dodge", stat="identity") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 60, hjust=1, size=13),
        axis.text.y = element_text(size=14),
        axis.title.x = element_text(size=14),
        axis.title.y = element_text(size=14),
        plot.title = element_text(size=14, hjust=0.5, face="bold"), legend.position = "none") +
  ggtitle("UK") +
  ylim(0,5) +
  xlab("") +
  ylab("Mean Support") +
  scale_x_discrete(labels= c("Social support", "Basic Income", "Job Guarantee", "Unskilled", "Skilled",
                             "Trade", "Reskill", "Tax automation"))

## Canada
summ <- multi %>%
  dplyr::filter(country=="Canada") %>%
  dplyr::group_by(group) %>%
  dplyr::summarize(socialsupport=mean(policy_socialspend, na.rm=T),
                   basicincome=mean(policy_basicincome, na.rm=T),
                   jobguarantee= mean(policy_jobguarantee, na.rm=T),
                   unskilled=mean(policy_unskilled, na.rm=T),
                   skilled=mean(policy_skilled, na.rm=T),
                   trade=mean(policy_trade, na.rm=T),
                   reskill=mean(policy_reskill, na.rm=T),
                   tax=mean(policy_automationtax, na.rm=T))

summ_long <- tidyr::gather(summ, policy_option, support, socialsupport: tax, 
                           factor_key=T)

CAN <- ggplot(summ_long, 
              aes(x=factor(policy_option), y=support, fill=group)) + 
  geom_bar(position="dodge", stat="identity") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 60, hjust=1, size=13),
        axis.text.y = element_text(size=14),
        axis.title.x = element_text(size=14),
        axis.title.y = element_text(size=14),
        plot.title = element_text(size=14, hjust=0.5, face="bold"), legend.position = "none") +
  ggtitle("Canada") +
  ylim(0,5) +
  xlab("") +
  ylab("Mean Support") +
  scale_x_discrete(labels= c("Social support", "Basic Income", "Job Guarantee", "Unskilled", "Skilled",
                             "Trade", "Reskill", "Tax automation"))

#figureA16
pdf(file = "figureA16.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
grid.arrange(AUS, CAN, UK, US, ncol=2)
dev.off()

## conjoint

#US

multi2_reshape_us <- dplyr::filter(multi2_reshape, country=="US")
#labels
attr(multi2_reshape_us$numhsafter, "label") <- "Number of High Skilled - After"
attr(multi2_reshape_us$numlsafter, "label") <- "Number of Low Skilled - After"
attr(multi2_reshape_us$wagelsafter, "label") <- "Wage of Low Skilled - After"
attr(multi2_reshape_us$wagehsafter, "label") <- "Wage of High Skilled - After"
attr(multi2_reshape_us$priceafter, "label") <- "Price - After"
attr(multi2_reshape_us$product, "label") <- "Product"
#marginal means for fairness

#figure A29
pdf(file = "figureA29.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(mm(data=multi2_reshape_us, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id), size=3) + theme(text = element_text(size = 15))   
dev.off()

# amce for fairness
# figure A30
pdf(file = "figureA30.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(cj(multi2_reshape_us, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id),
     size=3) + theme(text = element_text(size = 15))   
dev.off()

#marginal means for CEO
#figure A31
pdf(file = "figureA31.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(mm(multi2_reshape_us, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id), size=3) + theme(text = element_text(size = 15))  
dev.off()
# amce for CEO
pdf(file = "figureA32.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(cj(multi2_reshape_us, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id), size=3) + theme(text = element_text(size = 15)) 
dev.off()
#UK

multi2_reshape_uk <- dplyr::filter(multi2_reshape, country=="UK")
#labels
attr(multi2_reshape_uk$numhsafter, "label") <- "Number of High Skilled - After"
attr(multi2_reshape_uk$numlsafter, "label") <- "Number of Low Skilled - After"
attr(multi2_reshape_uk$wagelsafter, "label") <- "Wage of Low Skilled - After"
attr(multi2_reshape_uk$wagehsafter, "label") <- "Wage of High Skilled - After"
attr(multi2_reshape_uk$priceafter, "label") <- "Price - After"
attr(multi2_reshape_uk$product, "label") <- "Product"

#marginal means for fairness
#figure A25

pdf(file = "figureA25.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(mm(data=multi2_reshape_uk, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id), size=3) + theme(text = element_text(size = 15))   
dev.off()

# amce for fairness
#figure A26
pdf(file = "figureA26.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(cj(multi2_reshape_uk, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id),
     size=3) + theme(text = element_text(size = 15))   
dev.off()

#marginal means for CEO
#figure A27
pdf(file = "figureA27.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(mm(multi2_reshape_uk, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id), size=3) + theme(text = element_text(size = 15))  
dev.off()
# amce for CEO
# figure A28

pdf(file = "figureA28.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(cj(multi2_reshape_uk, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id), size=3) + theme(text = element_text(size = 15)) 
dev.off()

#Canada

multi2_reshape_can <- dplyr::filter(multi2_reshape, country=="Canada")
#labels
attr(multi2_reshape_can$numhsafter, "label") <- "Number of High Skilled - After"
attr(multi2_reshape_can$numlsafter, "label") <- "Number of Low Skilled - After"
attr(multi2_reshape_can$wagelsafter, "label") <- "Wage of Low Skilled - After"
attr(multi2_reshape_can$wagehsafter, "label") <- "Wage of High Skilled - After"
attr(multi2_reshape_can$priceafter, "label") <- "Price - After"
attr(multi2_reshape_can$product, "label") <- "Product"

#marginal means for fairness
#figure A21
pdf(file = "figureA21.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(mm(data=multi2_reshape_can, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id), size=3) + theme(text = element_text(size = 15))   
dev.off()

# amce for fairness
# figure A22
pdf(file = "figureA22.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(cj(multi2_reshape_can, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id),
     size=3) + theme(text = element_text(size = 15))   
dev.off()

#marginal means for CEO
#figure A23
pdf(file = "figureA23.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(mm(multi2_reshape_can, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id), size=3) + theme(text = element_text(size = 15))  
dev.off()

# amce for CEO
#figure A24
pdf(file = "figureA24.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(cj(multi2_reshape_can, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id), size=3) + theme(text = element_text(size = 15)) 
dev.off()

#Australia

multi2_reshape_aus <- dplyr::filter(multi2_reshape, country=="Australia")
#marginal means for fairness
attr(multi2_reshape_aus$numhsafter, "label") <- "Number of High Skilled - After"
attr(multi2_reshape_aus$numlsafter, "label") <- "Number of Low Skilled - After"
attr(multi2_reshape_aus$wagelsafter, "label") <- "Wage of Low Skilled - After"
attr(multi2_reshape_aus$wagehsafter, "label") <- "Wage of High Skilled - After"
attr(multi2_reshape_aus$priceafter, "label") <- "Price - After"
attr(multi2_reshape_aus$product, "label") <- "Product"

#marginal means for fairness
# figure A17
pdf(file = "figureA17.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(mm(data=multi2_reshape_aus, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id), size=3) + theme(text = element_text(size = 15))   
dev.off()

# amce for fairness
# figure A18
pdf(file = "figureA18.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(cj(multi2_reshape_aus, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id),
     size=3) + theme(text = element_text(size = 15))   
dev.off()

#marginal means for CEO
# figure A19
pdf(file = "figureA19.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(mm(multi2_reshape_aus, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id), size=3) + theme(text = element_text(size = 15))  
dev.off()
# amce for CEO
# figure A20
pdf(file = "figureA20.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(cj(multi2_reshape_aus, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
          priceafter + product, id= ~ response_id), size=3) + theme(text = element_text(size = 15)) 
dev.off()

### by group and by country

# fairness

# objective knowledge

#US

# amce


amces <- cj(multi2_reshape_us, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product + college, id= ~ response_id, estimate = "amce", 
            by = ~scenknow_all)

# diff in amce
diff_amces <- cj(multi2_reshape_us, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product + college, id= ~ response_id,estimate = "amce_diff", 
                 by = ~scenknow_all)
# figure A45
pdf(file = "figureA45.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()

#UK

#amce
amces <- cj(multi2_reshape_uk, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product + college, id= ~ response_id, estimate = "amce", 
            by = ~scenknow_all)
#diff in amce
diff_amces <- cj(multi2_reshape_uk, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product + college, id= ~ response_id,estimate = "amce_diff", 
                 by = ~scenknow_all)

# figure A41
pdf(file = "figureA41.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()

#Canada
#amce
amces <- cj(multi2_reshape_can, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product + college, id= ~ response_id, estimate = "amce", 
            by = ~scenknow_all)
#diff amce
diff_amces <- cj(multi2_reshape_can, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product + college, id= ~ response_id,estimate = "amce_diff", 
                 by = ~scenknow_all)
# figure A37
pdf(file = "figureA37.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()

#Australia
#amce
amces <- cj(multi2_reshape_aus, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product + college, id= ~ response_id, estimate = "amce", 
            by = ~scenknow_all)
#diff amce
diff_amces <- cj(multi2_reshape_aus, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product + college, id= ~ response_id,estimate = "amce_diff", 
                 by = ~scenknow_all)

# figure A33
pdf(file = "figureA33.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()


#self-knowledge

#US
#amce
amces <- cj(multi2_reshape_us, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product, id= ~ response_id, estimate = "amce", 
            by = ~self_know_automation_cat)
#diff amce
diff_amces <- cj(multi2_reshape_us, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product, id= ~ response_id,estimate = "amce_diff", 
                 by = ~self_know_automation_cat)

# figure A47
pdf(file = "figureA47.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()

#UK
#amce
amces <- cj(multi2_reshape_uk, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product, id= ~ response_id, estimate = "amce", 
            by = ~self_know_automation_cat)
#diff amce
diff_amces <- cj(multi2_reshape_uk, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product, id= ~ response_id,estimate = "amce_diff", 
                 by = ~self_know_automation_cat)
# figure A43
pdf(file = "figureA43.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()

#Canada
#amce
amces <- cj(multi2_reshape_can, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product, id= ~ response_id, estimate = "amce", 
            by = ~self_know_automation_cat)
#diff amce
diff_amces <- cj(multi2_reshape_can, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product, id= ~ response_id,estimate = "amce_diff", 
                 by = ~self_know_automation_cat)

# figure A39
pdf(file = "figureA39.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()

#Australia
#amce
amces <- cj(multi2_reshape_aus, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product, id= ~ response_id, estimate = "amce", 
            by = ~self_know_automation_cat)
#diff amce
diff_amces <- cj(multi2_reshape_aus, scenagree1 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product, id= ~ response_id,estimate = "amce_diff", 
                 by = ~self_know_automation_cat)

# figure A35
pdf(file = "figureA35.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()

# CEO

# knowledge

#US
#amce
amces <- cj(multi2_reshape_us, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product + college, id= ~ response_id, estimate = "amce", 
            by = ~scenknow_all)
# diff amce
diff_amces <- cj(multi2_reshape_us, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product + college, id= ~ response_id,estimate = "amce_diff", 
                 by = ~scenknow_all)

# figure A46
pdf(file = "figureA46.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()

# UK 
#amce 
amces <- cj(multi2_reshape_uk, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product + college, id= ~ response_id, estimate = "amce", 
            by = ~scenknow_all)
#diff amce
diff_amces <- cj(multi2_reshape_uk, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product + college, id= ~ response_id,estimate = "amce_diff", 
                 by = ~scenknow_all)

# figure A42
pdf(file = "figureA42.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()


# Canada
# amce
amces <- cj(multi2_reshape_can, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product + college, id= ~ response_id, estimate = "amce", 
            by = ~scenknow_all)
# diff amce
diff_amces <- cj(multi2_reshape_can, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product + college, id= ~ response_id,estimate = "amce_diff", 
                 by = ~scenknow_all)

# figure A38
pdf(file = "figureA38.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()

# Australia
# amce
amces <- cj(multi2_reshape_aus, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product + college, id= ~ response_id, estimate = "amce", 
            by = ~scenknow_all)
# diff amce
diff_amces <- cj(multi2_reshape_aus, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product + college, id= ~ response_id,estimate = "amce_diff", 
                 by = ~scenknow_all)

# figure A34
pdf(file = "figureA34.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()

#self-knowledge

#US
#amce

amces <- cj(multi2_reshape_us, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product, id= ~ response_id, estimate = "amce", 
            by = ~self_know_automation_cat)
#diff amce
diff_amces <- cj(multi2_reshape_us, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product, id= ~ response_id,estimate = "amce_diff", 
                 by = ~self_know_automation_cat)

# figure A48
pdf(file = "figureA48.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()

# UK
# amce
amces <- cj(multi2_reshape_uk, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product, id= ~ response_id, estimate = "amce", 
            by = ~self_know_automation_cat)
# diff amce
diff_amces <- cj(multi2_reshape_uk, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product, id= ~ response_id,estimate = "amce_diff", 
                 by = ~self_know_automation_cat)

# figure A44
pdf(file = "figureA44.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()

# Canada
# amce
amces <- cj(multi2_reshape_can, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product, id= ~ response_id, estimate = "amce", 
            by = ~self_know_automation_cat)

# diff amce
diff_amces <- cj(multi2_reshape_can, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product, id= ~ response_id,estimate = "amce_diff", 
                 by = ~self_know_automation_cat)

# figure A40
pdf(file = "figureA40.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()

# Australia
# amce
amces <- cj(multi2_reshape_aus, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
              priceafter + product, id= ~ response_id, estimate = "amce", 
            by = ~self_know_automation_cat)
# diff amce
diff_amces <- cj(multi2_reshape_aus, scenagree2 ~ numhsafter + numlsafter + wagehsafter + wagelsafter +
                   priceafter + product, id= ~ response_id,estimate = "amce_diff", 
                 by = ~self_know_automation_cat)

# figure A36
pdf(file = "figureA36.pdf",   # The directory you want to save the file in
    height = 8.27, # The width of the plot in inches
    width = 11.69) # The height of the plot in inches
plot(rbind(amces, diff_amces), size=3) + ggplot2::facet_wrap(~BY, ncol = 3L) + theme(text = element_text(size = 15)) 
dev.off()

## balance tables

# employed_short
multi <-
  multi %>%
  dplyr::mutate(employed = case_when(
    employment_screener=="Caring for family and working for pay" ~ "Employed",
    employment_screener=="Working for pay full-time" ~ "Employed",
    employment_screener=="Student and working for pay" ~ "Employed",
    employment_screener=="Working for pay part-time" ~ "Employed",
    employment_screener=="Retired" ~ "Retired",
    employment_screener=="Retired and working for pay" ~ "Retired",
    employment_screener=="Student" ~ "Student",
    employment_screener=="Self-employed" ~ "Employed",
    employment_screener=="Unemployed/looking for work" ~ "Unemployed",
    employment_screener=="Caring for family" ~ "Other",
    TRUE ~ "Other"
  ))

covs <- subset(multi, select = c(group, female,age, college, employed, income))
covs <- drop_na(covs, group,female)
covs$female <- ifelse(covs$female==1, "Female","Male")
covs$college <- ifelse(covs$college==1, "High Education","Low Education")
#balance <- bal.tab(covs, treat=covs$group)

tab_gender <- covs %>%
  group_by(group) %>%
  dplyr::count(female) %>%
  mutate(prop = prop.table(n))

tab_age <- covs %>%
  group_by(group) %>%
  summarize(age_mean=mean(age, na.rm=T),
            age_sd=sd(age, na.rm=T))


tab_college <- covs %>%
  group_by(group) %>%
  dplyr::count(college) %>%
  mutate(prop = prop.table(n))


tab_employed <- covs %>%
  group_by(group) %>%
  dplyr::count(employed) %>%
  mutate(prop = prop.table(n))

tab_income <- covs %>%
  group_by(group) %>%
  summarize(income_mean=mean(income, na.rm=T),
            income_sd=sd(income, na.rm=T))

pwc_income <- covs %>%
  pairwise_t_test(income ~ group, p.adjust.method = "bonferroni")
pwc_income
pwc_age <- covs %>%
  pairwise_t_test(age ~ group, p.adjust.method = "bonferroni")
pwc_age

pwc_employed <- chisq.test(covs$employed, covs$group)
pwc_female <- chisq.test(covs$female, covs$group)
pwc_college <- chisq.test(covs$college, covs$group)


cat <- dplyr::select(covs, group, female, college, employed)
cont <- dplyr::select(covs, group, age, income)

tbl_cat <- tbl_summary(cat, 
                       by = group, 
                       missing = "no",
                       statistic = all_categorical() ~ "{n} ({p}%)")

# Add chi-squared test results for categorical variables
tbl_cat <- add_p(tbl_cat, 
                 by = treatment, 
                 test = all_categorical() ~ "chisq.test", 
                 detailed = TRUE)

# Create a tbl_summary object for continuous variables
tbl_cont <- tbl_summary(cont, 
                        by = group, 
                        missing = "no",
                        statistic = all_continuous() ~ "{mean} ({sd})")

# Perform t-tests with Bonferroni correction for continuous variables
tbl_cont <- add_p(tbl_cont, 
                  by = group,
                  test = all_continuous() ~ "oneway.test",
                  detailed = TRUE)

tbl <- tbl_merge(list(tbl_cat, tbl_cont))

# Specify the columns to display in the table

# Convert the table to a 'gt' object
gt_tbl_cont <- as_gt(tbl_cont)
gt_tbl_cat <- as_gt(tbl_cat)

# Save the table to a file (e.g., html format)
gtsave(gt_tbl_cont, file = "table_cont.html")
gtsave(gt_tbl_cat, file = "table_cat.html")

